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The  pseudospectral  method  has  under-used  advantages  in  problems  involving  shocks  and 
discontinuities.  These  emerge  from  superior  accuracy  in  phase  and  group  velocities  as  com¬ 
pared  to  finite  difference  schemes  of  all  orders.  Dispersion  curves  for  finite  difference  schemes 
suggest  that  group  velocity  enor  typically  outranks  Gibbs'  error  as  a  cause  of  numerical 
oscillation.  A  flux  conservative  form  of  the  pseudospectral  method  is  derived  for  compatibility 
with  flux  limiters  used  to  preserve  monotonicity.  The  resulting  scheme  gives  high  quality 
results  in  linear  advection  and  shock  formation/propagation  examples.  ©  1989  Academic  Press. 
Inc. 


Introduction 


Spectral  methods  for  inviscid  flow  problems  involving  shocks  have  been 
developed  at  a  slower  rate  than  for  problems  where  discontinuities  are  absent.  This 
has  been  recognized  explicitly  [1]  and  implicitly  [2]  in  papers  taking  a  broad  view 
of  the  subject.  It  is  the  purpose  of  the  present  paper  to  show  how  the  pseudo¬ 
spectral  method  [3]  can  be  substituted  for  a  finite  difference  scheme  in  flux 
corrected  transport  (FCT)  calculations,  and  to  show  why  the  results  are  better.  The 
FCT  method  [4,  5]  is  effectively  a  local,  nonlinear  filter  used  to  control  numerical 
oscillations  near  shocks  and  steep  gradients. 

The  flux  correction  method  used  here  employs  the  Zalesak  flux  limiter  [5].  In 
this  approach,  two  arbitrary  schemes  are  cast  in  conservation  form:  one  of  high 
order  (pseudospectral),  and  the  other  monotone  (e.g.,  first  order  upwind).  A  “flux 
correction”  filter  effectively  assigns  a  local  weight  factor  to  fluxes  from  each  scheme. 
The  method  adopts  the  form  of  the  high  order  scheme  in  smooth  regions.  In  loca¬ 
tions  where  monotonicity  could  be  violated  by  the  high  order  scheme,  weighting  is 
shifted  in  favor  of  the  low  order  fluxes.  Comparisons  made  elsewhere  with  analytic 
solutions  show  that  the  accuracy  of  the  overall  method  generally  increases  with  the 
order  of  the  high  order  scheme  [6,  7]. 

Zalesak  carried  out  an  FCT  pseudospectral  calculatk^r  [6]  for  linear  advection 
which  gave  results  of  higher  quality  than  those  from  fyfite  difference  schemes  up  to 
sixteenth  order  in  space.  The  method,  however,  was  formulated  in  a  fashion  that 
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did  not  take  advantage  of  fast  Fourier  transforms  (FFTs).  Taylor  et  al.  [8]  gave 
results  of  a  pseudospectral  shock  calculation  using  FFT’s.  They  chose  a  Chebyshev 
polynomial  expansion  in  space,  appending  the  Boris-Book  flux  limiter  [4]  as  a 
postprocessor  to  control  oscillations  at  each  time  step.  This  involved  adding  global 
second  order  diffusion,  followed  by  flux  limited  antidifliision.  Properties  of  the 
method  itself  were  not  discussed  in  that  work,  nor  were  the  effects  of  wave  propaga¬ 
tion  through  the  nonuniform  grid  required  for  the  Chebyshev  expansion.  This  work 
takes  a  different  approach  to  discretization  and  flux  limiting  and  does  take  advan¬ 
tage  of  FFT’s. 

Depending  on  the  physical  problem  under  consideration,  high  order  methods 
may  have  advantages  over  shock  capturing  methods  of  the  Godunov  type  [9-1 1  ]. 
These  methods  discretize  variables  into  piecewise  polynomial  sections  (allowing 
discontinuities  between  sections  or  subsections)  and  use  gasdynamic  Riemann 
solutions  to  advance  a  step  in  time.  After  the  timestep,  the  solution  is  projected 
onto  the  original  basis  as  input  to  the  next  timestep.  These  methods  can  yield 
accurate  results  for  shocks,  but  their  complexity  increases  greatly  with  order.  In 
regions  where  the  flow  may  be  dominated  by  linear  advection,  these  methods  give 
results  equivalent  to  those  of  low  to  moderate  order  finite  differences. 


The  High  Order  Scheme:  Pseudospectral 

We  will  work  with  a  model  problem  for  a  scalar  hyperbolic  equation  in  one 
dimension: 

<5,d  +  <?x^(i;)  =  0.  (1) 

This  equation  leads  straight  to  the  conservation  of  J  v  dx  and  to  the  preservation  of 
monotonicity  (extrema  can  be  neither  created  nor  enhanced).  The  method 
developed  in  this  paper  may  be  modified  to  include  systems  of  equations  by 
replacing  <fi(v)  by  the  appropriate  advective  flux  uv,  and  using  suitable  donor  cell 
fluxes  in  Eq.  (10)  below. 

The  pseudospectral  method  evaluates  the  x  derivative  in  (1)  by  performing  an 
FFT,  then  evaluating  derivatives  analytically  before  the  inverse  FFT.  Let  the  finite 
difference  grid  be 

Xj=j-< 5x,  j=  1,2, ....  N,  (2) 

and  let  <f>  be  represented  in  space  by  Fourier  series: 

<0  =  L<A*)exp(ikx).  (3) 

* 

We  will  assume  periodic  boundary  conditions  on  all  variables  for  simplicity.  Other 
boundary  conditions  can,  of  course,  be  included  at  the  expense  of  altering  the 
selection  of  admissible  modes.  One  can  use  the  identity 


2 1 


(4) 
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to  express  the  x  derivative  of  (3)  as 


U  -ik{x  +  ix/2)  _  pik(x  -  Sxfl) 

SJ  =  Y,^(k)  Sin(/t  6x12)  '  (5) 

One  recognizes  this  as  the  difference  between  a  function  shifted  forward  half  a  grid 
cell  and  the  same  function  shifted  backward  half  a  grid  cell.  This  suggestive  result 
leads  directly  to  the  conservative  flux  form  given  in  (6)  and  (7)  below. 

For  the  discretization  of  ( 1 )  in  time,  we  will  use  the  leapfrog  trapezoidal  scheme 
[12]  with  step  St  between  time  levels.  This  scheme  can  be  thought  of  as  a 
predictor-corrector  with  the  predictor  being  a  leapfrog  step  from  time  level  n  —  1  to 
a  provisional  level  n  + 1  *  using  spatial  derivatives  evaluated  at  level  n.  Then  the 
provisional  level  n  +  1  *  is  averaged  with  n  to  yield  a  provisional  level  n  +  The 
trapezoidal  corrector  is  a  step  from  n  to  n  +  1  using  spatial  derivatives  evaluated  at 
n  +  j*.  This  provides  second  order  accuracy  in  time,  while  damping  the  computa¬ 
tional  mode  (“odd-even”  Nyquist  frequency  oscillations  in  space  and/or  time), 
which  could  interfere  with  the  flux  corrector  used  here. 

Employing  (5)  in  the  time-discretized  equation  (1)  leads  to  the  conservative  flux 
form  for  the  pseudospectral  method: 

«n+l  _ ..n  pn  +  1/2  ,  pn  +  1/2 

VJ  ~vj  j+i/2  +  fj-i/2’ 


where  superscripts  indicate  time  level  and 


pn  +  1/2 

V/  +  1/2 


=  —  Y  dn+ 1/2 (k)( emxl2  k6x/2  Y 
A  sin(&  Sx/2)J 


,ikx. 


(7) 


A  flux  form  of  this  type  seems  to  have  been  recognized  only  recently  as  a  possible 
ingredient  in  a  monotonicity  preserving  scheme  [13,  14],  The  leapfrog  trapezoidal 
time  advancement  scheme  used  here  is  stable  for  Courant  number 
e=  |u5t/5x|m„<21/2/7t  =  0.45,  where  u  is  the  characteristic  speed  (see  (9)  below). 
For  simple  leapfrog  advancement,  the  limit  is  l/n.  The  k  =  0  mode  has  been 
dropped  from  (7)  for  an  important  reason.  Different  additive  constants  in  high  and 
low  order  fluxes  would  lead  to  a  meaningless  bias  in  the  flux  correction  process, 
which  assigns  a  physical  interpretation  of  intergridpoint  transport  to  the  difference 
between  high  and  low  order  fluxes.  We  remove  the  spatial  mean  from  the  high 
order  flux  (7)  by  excluding  k  =  0,  and  from  the  low  order  flux  by  computing  the 
average  and  subtracting  it. 


The  Low  Order  Scheme:  Monotone  Upwind 

A  convenient  first  order  monotone  scheme  [7]  for  (1)  is  upwind  differencing  with 
the  differencing  direction  defined  by  a  variable  Tj 

WJ+  1/2  =  (<t>j+  I  —  4>j)  •  (Vj+  1  —  Vj),  (8)  - 

/ I  I  . 


on 
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whose  sign  is  that  of  the  characteristic  speed  u  somewhere  on  the  interval  (vjt  vJ+l), 
where 


u  — 


dj_ 

dv 


We  define  first  order  fluxes 

_j<pjSt/5x, 
Jj+U2  1  +J+,St/Sx, 

and  remove  the  spatial  average, 


wj+  1/2^0 
W'/+l/2<0, 


fj  +  1/2  ~*fj  +  1/2  (/)• 


(9) 

(10a) 


(10b) 


When  equations  are  given  without  superscripts,  it  is  understood  that  all  variables 
are  evaluated  at  the  same  time  level.  The  removal  of  the  mean  value  in  (10b)  has 
no  effect  on  the  low  order  update  (11),  but  is  necessary  for  compatibility  with  (7) 
in  the  flux  correction  process.  A  monotone  first  order  update  accounting  for 
reversals  in  the  sign  of  u  is 


»7+,'-»;-/;+12+/7-«2-  (in 

This  provisional  update  to  time  level  n  4-  1'  is  set  aside  for  correction  as  described 
below.  The  above  scheme  is  monotone  for  e  ^  j.  This  allowed  Courant  number 
is  smaller  by  a  factor  of  two  than  that  for  stability  alone.  The  extra  factor  of 
two  allows  for  monotonicity  preservation  when  u  is  compressive  at  an  extremum 
of  v  [7], 


Flux  Correction:  A  Nonlinear  Filter 


The  next  step  is  to  construct  a  set  of  “flux  corrections” 


S/»  +12 _  £vi  +  I  2  fn 

j+  1/2  —  r  j+  12  i  +  i/2 


(12) 


and  filter  them  with  a  flux  limiter  to  prevent  creation  or  enhancement  of  extrema. 
For  problems  of  higher  than  one  dimension,  flux  limiters  have  been  constructed  [5] 
which  require  auxiliary  storage.  In  one  dimension,  however,  it  is  sufficient  to  use 


+  1/2 
°J  j+  1/2 


•  Sj+  u2  ■  max{0,  min[7}+  Sj  + I/2  •  K+V  -  v”*tv). 


sJ+l/2-(vri-vp!m 


(13) 


where 


■5/+ 1/2  —  sign(<5/7  *,'£), 

2y+ 1/2  =  1/2  4"  sign(y"+ ('  —  y"+  1  )]. 


;  1 

i 
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This  is  similar  to  the  Boris-Book  flux  limiter,  but  with  the  following  differences.  The 
flux  correction  Sf  has  been  constructed  as  the  difference  between  fluxes  from  two 
specified  schemes,  rather  than  the  negative  of  a  diffusive  flux  added  to  damp  oscilla¬ 
tions.  The  factor  TJ+ 1/2<5/"+ 1V2  replaces  I  <5/7+ I  'n  tlie  limiter  of  Ref.  [1].  This 
replacement  requires  that  the  correction  to  high  order  not  decrease  gradients 
locally;  i.e.,  the  flux  correction  should  be  of  an  antidiffusive  nature  everywhere. 

This  form  preserves  monotonicity  of  v  and  results  in  total  fluxes  which  are  point- 
by-point  intermediate  between  high  and  low  order  fluxes.  The  completed  monotone 
pseudospectral  update  is  then 

vTj  +  '  -»;  +  r -SfX$  +  6ffr,g.  (15) 

This  approach  to  flux  limiting  recovers  the  high  order  fluxes  [5]  FJ+ 1/2  at  all 
locations  j+  1/2  except  where  the  high  order  scheme  could  violate  monotonicity  at 
j  or  7+  1.  Thus  spectral  accuracy  is  retained  in  smooth  regions. 


Oscillations:  Gibbs'  Error  and  Phase  Error 

This  section  discusses  properties  of  numerical  schemes  which  have  been 
recognized  previously  [15-17],  but  which  have  received  little  attention  in  the  con¬ 
struction  of  flux  limiters.  Numerical  oscillation  is  commonly  attributed  to  Gibbs’ 
error;  i.e.,  the  tendency  of  a  series  representation  to  oscillate  near  a  jump  in  nodal 
values.  A  stronger  contributor  to  numerical  oscillation  in  advection  scheme  results 
is  discretization  error  in  phase  and  group  speed.  Numerical  oscillations  from 
difference  schemes  of  various  orders  are  illustrated  in  Fig.  1.  (This  figure  is  the  first 
part  of  a  before-and-after  pair,  with  Fig.  4  containing  results  after  flux  correction. 
The  pseudospectral  results  in  Fig.  4  are  flux-corrected  in  exactly  the  same  way  as 
the  finite  difference  results.)  Linear  advection  examples  are  given  for  two  profiles 
which  reveal  numerical  artifacts  in  finite  difference  schemes:  a  discontinuous  square 
wave  and  a  semicircle.  The  square  wave  tests  monotonicity  preservation  and  reveals 
residual  numerical  diffusion  in  the  rounding  of  the  discontinuity.  The  semicircle  is 
sensitive  to  terracing  or  oscillation  in  slope. 

A  trademark  of  Gibbs’  error  is  that  the  width,  rather  than  the  amplitude  of  the 
oscillation,  tends  to  zero  as  the  number  of  retained  frequencies  increases.  The  width 
of  the  region  of  spurious  oscillations  in  Fig.  1,  however,  increases  with  the  order  of 
the  scheme.  This  is  characteristic  of  group  velocity  error.  Figure  2  gives  a  result 
more  typical  of  Gibbs’  error.  Here  the  linear  advection  example  of  Fig.  le  is  carried 
much  further  in  time  with  the  pseudospectral  method.  A  mild  damping  of  the  com¬ 
putational  mode  results  from  the  leapfrog-trapezoidal  time  advancement  scheme.  As 
pointed  out  below,  the  pseudospectral  method  gives  accurate  phase  speeds  for  all 
modes  except  the  odd-even  mode.  Gradual  removal  of  the  small  spectral  interval 
where  phase  and  group  velocity  errors  are  significant  has  exposed  the  underlying 
Gibbs’  error. 


Fig.  1.  Passive  scalar  advection  test  without  flux  limiter,  (a)  Initial  condition;  (b)  Second  order 
spatial  differencing  after  30  timesteps  at  a  Courant  number  of  0.1;  (c)  Fourth  order;  (d)  Eighth  order; 
(e)  Pseudospectral.  Solid  lines:  analytic  solution;  points:  numerical  results. 
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Figure  1 — Continued 

For  the  sake  of  illustration,  Fig.  3  shows  phase  and  group  velocities  for  eight 
schemes  of  order  2  through  128  plus  pseudospectral.  Such  high  orders  are  included 
only  to  illustrate  how  much  would  be  required  of  a  finite  difference  scheme  to 
challenge  pseudospectral  accuracy  at  high  wavenumber.  These  curves  are  calculated 
to  arbitrary  order  as  follows.  Consider  the  linear  advection  equation 

d,v=—cdxv,  (16) 

where  c  is  constant.  Assuming  temporal  and  spatial  variation  ei{kx~w,)  and 
replacing  dx  in  (16)  with  a  finite  diffeence  operator  of  a  given  order  accuracy  then 
gives  co  as  a  function  of  k  dx.  (Temporal  truncation  error  is  omitted  from  this 
analysis.)  For  a  stable  scheme  of  even  order,  co  is  real.  An  expression  for  the 
numerical  derivative  of  a  function  f(x)  accurate  to  order  2 M  at  a  point  x0  on  an 
evenly  spaced  grid  is  [18] 

8f(x0)  y  f(x,)-f(x  o)  j 

8x  +  0  Xi~X0  j-  -M  +  0,i  j~i 


(17) 
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Fig.  4.  Results  analogous  to  Fig.  1  after  800  steps  including  the  flux  limiter  of  Eq.  (13). 
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Figure  4 — Continued 


Figure  3  reflects  that  phase  velocities  from  finite  difference  schemes  are  accurate 
for  the  well-resolved  (low  wavenumber)  modes,  but  fall  to  zero  at  the  odd-even 
mode.  Then  the  best  a  real  scheme  can  do  is  stay  with  the  linear  dispersion  relation 
as  long  as  possible  before  falling  away.  (In  principle,  a  complex  scheme  could  be 
constructed  to  retain  phase  information  lost  in  time-advancing  the  odd-even 
mode.)  Since  the  group  velocity  (i.e.,  the  wave  packet  speed)  is  proportional  to  the 
slope  of  the  dispersion  curve  co(fc),  information  in  the  high  wavenumber  modes 
propagates  in  the  wrong  direction  at  a  speed  which  increases  with  increasing  order. 
This  behavior  is  evident  in  Fig.  1. 

With  the  pseudospectral  method,  only  the  highest  mode  falls  off  the  phase 
velocity  curve,  so  that  group  velocities  are  accurate  for  all  but  the  highest  mode 
interval.  The  leapfrog-trapezoidal  time  advancement  scheme  used  here  contributes 
a  mild  damping  to  these  modes,  keeping  results  from  being  highly  contaminated 
with  the  backward-running  error  seen  in  Fig.  1.  Extension  of  the  calculation  of 
Fig.  le  from  30  to  800  timesteps  with  the  pseudospectral  method  (Fig.  2)  yields  the 
low  oscillation  and  high  structure  resolution  levels  one  would  hope  for  before 
applying  a  monotonicity  filter. 

An  interesting  and  little-realized  fact  can  be  gleaned  from  Fig.  3.  Even  when  the 
finite  difference  order  is  one  less  than  the  number  of  grid  points,  errors  still  occur 
in  the  high  wavenumber  phase  and  group  velocities,  and  the  pseudospectral  method 
stands  unsurpassed.  In  this  case,  the  polynomial  representing  the  finite  difference 
scheme  passes  through  all  the  nodal  values,  and  one  might  be  tempted  to  think  that 
this  situation  could  not  be  improved.  The  resolution  of  this  apparent  paradox  is 
that  finite  differences  do  not  give  exact  results  for  differentiation  of  the  sinusoid 
used  in  defining  phase  velocity,  while  the  pseudospectral  method  does. 
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Figure  5 — Continued 


Results  with  Flux  Correction 


Linear  Advection 

Figure  4  shows  results  analogous  to  those  of  Fig.  1  with  a  flux  limiter  of  the  form 
of  ( 1 3 )  applied  to  finite  difference  methods  of  increasing  order,  including  the  scheme 
of  (7 )— ( 13).  This  latter  combination  will  be  designated  PSF  for  pseudospectral  flux 
correction.  All  tests  use  the  leapfrog-trapezoidal  time  differencing  scheme.  Lacking 
the  spurious  oscillations  of  Fig.  2,  Fig.  4  shows  the  kind  of  convergence  with 
increasing  order  that  one  would  hope  to  find  in  a  well-behaved  scheme.  Numerical 
smoothing  evident  in  the  rounding  and  broadening  of  the  step  function  profile 
decreases  with  order,  even  after  the  addition  of  the  local  diffusion  implicit  in  the 
flux  limiter. 

A  more  demanding  test  of  the  schemes  than  monotonicity  preservation  is  their 
accuracy  against  the  formation  of  terraces.  These  result  from  the  tendency  to  form 
oscillations  in  slope.  These  slope  oscillations  could  occur  at  low  amplitude  without 
invoking  any  local  smoothing  from  the  flux  limiter.  Only  when  they  reach  an 
amplitude  large  enough  to  cause  artificial  extrema,  the  flux  limiter  enters  and 
effectively  converts  an  individual  artificial  extremum  to  a  plateau.  A  distribution  of 
these  plateaus  appears  as  a  terrace.  The  semicircle  profile  of  Fig.  4  shows  a 
decreasing  tendency  to  form  terraces  as  the  order  of  the  scheme  increases.  Both  the 
amplitude  and  width  of  the  terraces  decrease  with  order.  This  is  consistent  with  the 
view  that  increasing  the  order  of  the  scheme  moves  to  a  higher  wavenumber  and 
narrows  the  spectral  interval  over  which  phase  and  group  velocities  are  in  error. 

Shock  Formation 

Another  revealing  test  of  a  scheme  for  hyperbolic  problems  is  to  follow  the  evolu¬ 
tion  of  a  continuous  profile  subject  to  the  inviscid  Burgers’  equation  into  steepening 
and  shock  propagation.  Many  tests  of  contemporary  shock  following  methods  have 
been  published,  in  which  shocks  are  present  in  initial  conditions.  These  tests  may 


426 


FLUX-CORRECTED  PSEUDOSPECTRAL  METHOD 


427 


Figure  6 — Continued 


not  reveal  artifacts  which  can  occur  in  the  continuously  steepening  case. 
Pathologies  that  can  occur  in  steepening  have  been  demonstrated  elsewhere  [7]. 

Figure  5  shows  PSF  results  for  shock  formation  and  propagation  as  described  by 
the  inviscid  Burgers’  equation,  in  comparison  with  the  analytic  solution.  The  shock 
front  behavior  is  in  excellent  agreement  with  the  analytic  solution,  while  the  back 
side  of  the  shock  is  contaminated  by  an  odd-even  oscillation.  This  is  substantially 
improved  in  Fig.  6  by  the  use  of  a  high  order  (16th)  linear  filter.  The  filter  consists 
of  adding  to  (7)  a  flux 

DJ+}/}  =  -eAm/2  -  >;+ ,  -  v"),  (20) 

where  e  is  the  maximum  Courant  number  in  the  calculation,  m=  16,  and  A  is  the 
normalized  second  difference  operator, 

AVj 3  -  ivJ+ ,  +  {Vj -  Jt>,._  , .  (21 ) 

The  filter  (20)  is  added  to  the  high  order  flux  (7)  before  the  flux  limiter  is  evoked, 
so  that  monotonicity  of  the  total  scheme  is  preserved.  This  form  for  the  filter  was 
found  most  satisfactory  after  experimenting  with  several  more  complicated  forms. 


Summary 

Gibbs’  error  is  often  cited  as  the  cause  of  numerical  oscillation  near  discon¬ 
tinuities.  It  is  likely,  however,  that  finite  difference  phase  error  causes  more  artificial 
oscillation  than  Gibbs’  error.  The  pseudospectral  method  is  capable  of  higher 
quality  resolution  of  discontinuities  than  finite  difference  schemes  because  of  its 
accuracy  in  phase  and  group  velocities  for  spatial  modes  below  the  Nyquist 
frequency.  Numerical  oscillations  in  both  value  and  slope  are  smaller  with  the 
pseudospectral  method  than  with  finite  difference  schemes.  The  pseudospectral 
method  when  cast  in  flux  difference  form  (7)  is  compatible  with  monotonicity- 
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enforcing  flux  limiters  developed  for  finite  difference  schemes.  The  resulting 
pseudospectral  flux  correction  (PSF)  method  appears  very  attractive  for  problems 
involving  discontinuities  and  shock  formation.  For  passive  advection  of  a 
discontinuity  at  a  uniform  speed,  PSF  controls  oscillations  in  value  and  slope  while 
maintaining  structure  details.  For  shock  propagation,  however,  a  high-order  linear 
filter  is  desirable  for  controlling  oscillations  in  slope  behind  the  shock.  Further 
work  is  needed  on  this  feature  of  the  scheme. 
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